function [inequ_par, tau_hsv, lambda_hsv] = calibrate_inequity_aversion(parameter_candidates, w0, wf, knitro_path)

    addpath(knitro_path);

    addpath('../')
    addpath('../Densities')
    addpath('../UtilityFunctions')
    addpath('../ProductionFunctions')
    addpath('../Economies/')

    calibration_version = '20191107_nochoice';
    prefix = calibration_version;
    prefix_pf = prefix;
    prefix_skill = prefix;

    dir_in = '../../../data/generated/matlab/Calibration/';
    dir_out = dir_in;

    const = tools.csv2struct('../../../data/raw/constants.csv');

    % load calibrated parameters for skill distribution
    par_dist = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_sol.csv']);
    mom_dist = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_mom.csv']);

    DensSkill = abilityDistFixed(par_dist);
    DensPart = ParticipationNormal(par_dist);

    % utility function object
    par_util.epsilon = const.epsilon;
    Util = UtilityQuasiLinear(par_util);

    % combine in SkillAndParticipation object
    par_dist.dens_skill = DensSkill;
    par_dist.dens_part = DensPart;
    par_dist.util = Util;

    par_dist.transfer = par_dist.transfer;
    par_dist.num_nodes_phi = 100;

    Dist = SkillAndParticipation(par_dist);

    % Production function 
    par_prod_fun =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_sol.csv']);
    mom_prod_fun =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_mom.csv']);

    par_prod_fun.parse_inputs = true;
    par_prod_fun.KrusellRobotsEach_use_properties = true;
    ProdFun = KrusellRobotsEach(par_prod_fun);

    % Economy object
    par_econ.varnames = {'L_M','L_R','L_C', ...
                        'K_B_M','K_B_R','K_B_C', 'K_E','K_S',...
                        'tau_B','tau_E','tau_S','transfer','tau_hsv','lambda_hsv','compensation'};
    par_econ.pf = ProdFun;
    par_econ.util = Util;

    par_econ.lb = 1e-16 .* ones(1,8);
    par_econ.ub = 1e16 .* ones(1,8);

    % robot prices and taxes
    inp.tau_B = par_prod_fun.tau_B;
    inp.tau_E = par_prod_fun.tau_E;
    inp.tau_S = par_prod_fun.tau_S;

    par_econ.q_B = mom_prod_fun.q_B;
    par_econ.q_E = mom_prod_fun.q_E;
    par_econ.q_S = mom_prod_fun.q_S; 

    q_B_0 = par_econ.q_B;

    par_econ.dist = Dist;

    par_econ.solver = 'knitro_eqsolve';

    par_econ.phi_lb = const.phi_lb;
    par_econ.par_dist = par_dist;

    par_econ.w_lb = w0;
    par_econ.w_ub = wf;

    econ = KrusellWithRobotsEachEconomy(par_econ);
    econ.set_num_nodes(100,50,20,50);

    % compute mass for normalization
    integrand = @(w) econ.dist.dens_skill.f(w,par_dist.Y_M_0,par_dist.Y_R_0,par_dist.Y_C_0,par_dist);
    mass = integral(integrand,w0,wf);

    econ.dist.dens_skill.mass = mass;

    integrand = @(w) econ.dist.dens_skill.f(w,par_dist.Y_M_0,par_dist.Y_R_0,par_dist.Y_C_0,par_dist);
    mass_normalized = integral(integrand,w0,wf)

    inp.L_M = par_prod_fun.L_M_0 + 1e-3;
    inp.L_R = par_prod_fun.L_R_0 + 1e-3;
    inp.L_C = par_prod_fun.L_C_0 + 1e-3;
    inp.K_B_M = par_prod_fun.K_B_M_0 + 1e-3;
    inp.K_B_R = par_prod_fun.K_B_R_0 + 1e-3;
    inp.K_B_C = par_prod_fun.K_B_C_0 + 1e-3;
    inp.K_E = par_prod_fun.K_E_0;
    inp.K_S = par_prod_fun.K_S_0;
    inp.transfer = par_dist.transfer;
    inp.tau_hsv = par_dist.tau_hsv;
    inp.lambda_hsv = par_dist.lambda_hsv;
    inp.compensation = 0;

    econ.revenue_requirement = mom_prod_fun.revenue_requirement;

    lb.L_M = 1e-6;
    lb.L_R = 1e-6;
    lb.L_C = 1e-6;
    lb.K_B_M = 1e-12;
    lb.K_B_R = 1e-12;
    lb.K_B_C = 1e-12;
    lb.K_E = 1e-6;
    lb.K_S = 1e-6;
    lb.tau_B = inp.tau_B;
    lb.tau_E = inp.tau_E;
    lb.tau_S = inp.tau_S;
    lb.transfer = 0;
    lb.tau_hsv = 0.05;
    lb.lambda_hsv = 1;

    ub.L_M = 1e6;
    ub.L_R = 1e6;
    ub.L_C = 1e6;
    ub.K_B_M = 1e8;
    ub.K_B_R = 1e8;
    ub.K_B_C = 1e8;
    ub.K_E = 1e6;
    ub.K_S = 1e6;
    ub.tau_B = inp.tau_B;
    ub.tau_E = inp.tau_E;
    ub.tau_S = inp.tau_S;
    ub.transfer = 100;
    ub.tau_hsv = 1;
    ub.lambda_hsv = 1e5;    
    
    lb.compensation = 0;
    ub.compensation = 0;
    
    tmp = arrayfun(@(x) objective(x), parameter_candidates, ...
                      'UniformOutput', false);
    
    mat_tmp = cell2mat(tmp);
    
    [~, col_idx] = min(mat_tmp(1,:));
    inequ_par = parameter_candidates(col_idx);
    tau_hsv = mat_tmp(2,col_idx);
    lambda_hsv = mat_tmp(3,col_idx);

    store.inequ_par = inequ_par;
    tools.struct2csv(store, [dir_out,'inequ_par.csv']);
    
    fprintf('\ncalibrated inequity aversion: %f at tau_hsv: %f, lambda_hsv: %f\n',inequ_par, tau_hsv, lambda_hsv)
     
    function out = objective(inequity_aversion_parameter)
        econ.r = inequity_aversion_parameter;
        [sol,fval,exitflag,output] = econ.compute_optimtax(inp,lb,ub);
        tmp = [sol.tau_hsv - inp.tau_hsv,...
              sol.lambda_hsv - inp.lambda_hsv];
        ssr = tmp * tmp';
        lambda_hsv = sol.lambda_hsv;
        tau_hsv = sol.tau_hsv;
        out = [ssr,tau_hsv,lambda_hsv]';
    end
        
end